A study of a WENO-TVD finite volume scheme for the numerical 
simulation of atmospheric advective and convective phenomena 



Dante Kalise^ 

^Dipartimento di Matematica, Universitd di Roma "La Sapienza", P. Aldo Mow 2, 00185 Roma, Italy. 



Abstract 

We present a WENO-TVD scheme for the simulation of atmospheric phenomena. The scheme 
considers a spatial discretization via a second-order TVD flux based upon a flux-centered limiter 
approach, which makes use of high-order accurate extrapolated values arising from a WENO 
reconstruction procedure. Time discretization is performed with a third order RK-TVD scheme, 
and splitting is used for the inclusion of source terms. We present a comprehensive performance 
study of the method in atmospheric applications involving advective and convective motion. We 
present a set of tests for space-dependent linear advection, where we assess convergence and 
robustness with respect to the parameters of the scheme. We apply the method to approximate 
the 2D Euler equations in a series of tests for atmospheric convection. 

Keywords: WENO reconstruction, TVD schemes, Runge-Kutta methods, splitting, limiter, 
centered schemes, swirling flow, frontogenesis, Euler equations, convection 



1. Introduction 

Advection and convection, understood as the class of phenomena related to fluid flow motion 
or transport is a concept covering difl'erent relevant situations in atmospheric modelling. The 
convention in the literature is to use the term advection when a quantity experiences motion due 
to the presence of an acting velocity field, which is most often related to horizontal motion, while 
convection refers to motion caused by thermodynamic considerations, primarily occurring in the 
vertical direction. In this article we are concerned with the development of an accurate scheme 
able to handle both types of motion. Models which describe such behavior are, in a first step, 
linear scalar advection models, including space-dependent velocity fields and, in a more elabo- 
rated formulation, the set of Euler equations for gas dynamics, parameterized in scales that are 
typical of atmospheric phenomena. Achieving an accurate and physically meaningful numerical 
approximation of such models is undoubtedly, a challenging task. During the last decades, a 
method which has gained popularity among the atmospheric modelling community is the finite 
volume (FV) framework [Jj i2j| . In a computationally efficient way, it preserves many aspects 
of the underlying physics (as conservation for instance), while allowing formulations that are 
able to achieve a high level of accuracy, which is fundamental for numerical weather prediction. 
In the FV context there is a considerable amount of available methods oriented to the solution 
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of hyperbolic system of conservation laws. We are concerned with a particular class of those 
methods, the so-called WENO (weighted essentially non-oscillatory methods), which is a class 
of methods allowing the generation of high-resolution approximations in space via a polyno- 
mial reconstruction procedure from cell averages | 3 , 4 |; this technique is combined with suitable 
fluxes and time marching procedures in order to generate a method of global high-order of accu- 
racy. Moreover, at every step, by adequately enforcing the concept of total variation diminishing 
(TVD), the method avoids the production of spurious oscillations and preserves monotonicity. 
This article addresses every step in the construction of the above described scheme. Once the 
WENO reconstruction of the variables is performed, high-order extrapolated values are avail- 
able for calculation of numerical fluxes across the cell interfaces. Based upon f5\, we make use 
of a second-order TVD flux based on a centered limiter (FLIC) approach [6 1 ; we highlight the 
use of the centered approach instead of classical upwinding considerations |7 |. After the spatial 
discretization routine is completed, the scheme moves forward in time via a Runge-Kutta TVD 
method |^ with ensures stability, preservation of the spatial accuracy and avoids the generation 
of new extremal values. An important part of this work is devoted to the numerical validation 
of the proposed scheme. We first address 2D models of space-dependent advection where sev- 
eral features of the scheme are tested, to then study convective phenomena based upon the Euler 
equations, with a set of well-known tests for atmospheric modelling. The performance of the 
scheme is assessed in terms of accuracy, its ability to preserve monotonicity in the presence of 
sharp solutions, its robustness with respect to flux parameters, correct front locations and energy 
conservation. 

The article is structured as follows. In section 2 we present the full numerical WENO-TVD 
scheme. In section 3 we develop a comprehensive study of the method for advective models, 
while section 4 is devoted to the analysis of convective phenomena. Final remarks are discussed 
in section 5. 



2. A numerical scheme for the system of balance laws 

In this section we present a finite volume scheme for a two-dimensional system of balance 
laws of the form 

dtQ^ d.r{Q)+ d,^{Q)= S{Q\ (1) 

where 2 is a vector of conserved variables, T and are physical fluxes, and <S is a source 
term. We first indicate that our strategy will be based in a splitting scheme, as it is suggested in 
191 given the flux choice that we will make. Thus, we will first establish a numerical scheme for 
the system of conservation laws 

dtQ+ d,r(Q)^ d,9{(Q) = 0, (2) 

to be combined with a procedure for the resolution of the source term dynamics 

dt(Q)=S(Q). (3) 

In order to approximate eq. ([2]) we begin by meshing the spatial domain flx,z ii^to uniform control 
volumes Q/j = l^i-i/i, ^i+i/i] x [Zj-i/2, Zj+1/2] of size AxAz; inside every control volume we 
average with respect to x and z leading to the semi-discrete scheme 

dQij(t) _ 1 1 

= -^(^/+l/2j - Fi-i/2,j) - —(Hij+i/2 - Hij-1/2) = Lij(Q)^ (4) 
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where 

I I rxi+m rzj+i/2 
Qi'j=irir Q{x,z,t)dzdx, (5) 

Fm/2j= ^ f' ' r{Q{xM,2,z,t))dz, Hij^i,2 = ^f ' 'H(Q{x,Zj+i/2,t))dx. (6) 
We approximate the expressions in (|6]( by conventional Gaussian quadrature formulas 

Pi+\l2,j ^ZGp ^ ZGp, 0) dz , Hij+ij2 ^Q^^Gp,Zj+l/2, 0) dz, 

ZGp XGp 

(7) 

where xop and zop are prescribed Gauss points with corresponding weights Wx^p and w^^p respec- 
tively. The computation of ^ is performed via a high-resolution approach that makes use of a 
WENO reconstruction procedure; after this step is completed, a polynomial of prescribed order 
is obtained at every cell, and therefore, at every cell interface, accurate flux calculations can be 
performed by taking extrapolated boundary values. 

We briefly describe the WENO reconstruction procedure that is used in this article; we opted 
for the technique described in |[lOl[TTl in its third order (quadratic reconstruction) version. This 
technique makes extensive use of the structure of the reconstruction procedure in one dimension, 
adding some additional mixed terms ("cross terms") that are efl&ciently computed by reduced 
stencils. It is an optimal and easy way to implement the algorithm for achieving high-order 
reconstructions in 2 and 3 dimensions; it also defines an unique polynomial in every cell, which 
is particularly useful when space dependent source terms such as viscosity are considered. At 
a given time t (the subscript indicating time is omitted throughout this derivation), given the 
set of averaged values [Qij] for the whole domain, at every cell, the reconstruction procedure 
seeks a quadratic expansion upon a linear combination of Legendre polynomials rescaled in 
local coordinates (x,z) = [-1/2, 1/2] x [-1/2, 1/2] expressed in the form 

e(x, z) = Go + QxPi{x) + QMx) + Q,Pi{z) + QzzPiiz) + QxzPi{x)Pi{z\ (8) 



Pi{x) = x P2(x) = x^-^. (9) 

Except for the last term in ([5]), every coefficient can be computed by performing a dimension-by- 
dimension reconstruction, which we now illustrate. We assign the subscript "0" to the cell where 
we are computing the coefficients, other values indicating location and direction with respect 
to Qo (note that the notation is coherent with the fact that the first coefficient in the expansion 
Qo, holds Qo = Qij, i.e., the centered value). Next, for this particular problem we define three 
stencils 

S' ={Q.2M-uQo}, S^ = {Q.uQoMi}, = {Qo,QuQ2}, (10) 
and in every stencil we compute a polynomial of the form 

Q^Hx) = Qf + QfPi(x) + QfxPiix) / = 1, 2, 3. (1 1) 
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The coefficients are given by 

S' : Q(^i> = -2Q_i + e_2/2 + 320/2, 2^ = (2-2 - 2e_i + eo)/2 , (12) 

: ef = (Qi - e-i)/2, = (2-1 -2Qo + 20/2, (13) 

■■ 2?' = -32o/2 + 22i-22/2, 2^:' = (Co - 22-1 + 22)/2. (14) 
For every polynomial we calculate a smoothness indicator defined as 

/5® = (2«f + (15) 
leading to the following WENO weights: 

^^'^=^ ' ^^'^= . \r.. ^ (16) 

where 6 is a parameter introduced in order to avoid division by zero; usually 6 = 10"^^. The 



scheme is rather insensitive to the parameter r, which we set r = 5. The parameter A is usually 
computed in an optimal way to increase the accuracy of the reconstruction at certain points; we 
opt for a centered approach instead, thus /l^^^ = /l^^^ = 1, while /l^^^ = 100. The ID reconstructed 
polynomial is given by 

Q(x) = oj^^^Q^^\x) + aj^^^Q^^\x) + aj^^^Q^^\x). (17) 

Next, a ID reconstruction in the z direction is performed in a totally analogous way. Finally, we 
address the computation of the mixed term Q^z, which is calculated in a 2D fashion. Keeping the 
same convention regarding location subscripts as in ID, 1 10] considers 4 formulas for the cross 
term upon taking all the moments around the cell. The expressions for the cross term are: 



2(1) = 


-- Qui - 2o,o -Qx-Qz- Qxx - Qzz, 


(18) 


2(2) . 


-- -2i,-i + 2o,o + 2x - 2<: + Qxx + Qzz' 


(19) 


Q^S - 


= -2-1,1 + 2o,o - 2x + 2<: + Qxx + Qzz' 


(20) 


Q'S - 


-- 2-1,-1 - Qo,o + Qx + Qz- Qxx - Qzz, 


(21) 



and the corresponding smoothness indicators are given by 

/5«=4(2«)%4(2«f + (2«)'. (22) 

Note that in the first part of the reconstruction, when the weights were computed, a larger subop- 
timal weight was assigned to the central stencil, which is a way to ensure stability and robustness 
of the algorithm by sacrificing additional order in the approximation (for more details, see CD). 
However, for this term, the numerators assigned to the corresponding a's remains the same for 
every expression. The computation of this term concludes the reconstruction procedure, and now 
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we have at our disposal one polynomial per cell that can be used to calculate values at the bound- 
aries or inside the cell. The next step in our numerical scheme consists of the calculation of the 
numerical fluxes (|7]), which will use extrapolated boundary values of the reconstructed polyno- 
mials. Rather than the use of the classical WENO scheme (as in |3 | or |4|), which performs 
this calculation via a first order flux, we opt for the WENO-TVD approach described in |5|. We 
make use of the 2D extension of the flux-limiter-centred scheme (FLIC) approach presented in 
161 UJ, which is a second-order, centered and non-oscillatory flux. In our case, it consists of a 
flux-limited version of a generalized Lax Wendrofl' flux, using as a low-order flux the GFORCE 
(generalized first order centred) flux |[T3l , which can be interpreted as a convex combination of 
Lax-Friedrichs and Lax-Wendrofl'-type of fluxes: 



where 



^GFORCE 
i+l/2,j 



T^FLIC _ T^GFORCE , 
^i+l/2,j - ^i+l/2,j ^ 



J. (t^lw t^GFORCE\ 



LF 

1+1/2 J ' 



(23) 



(24) 



i+U2,j 



(26) 



The parameter oj varies between and 1, and is chosen in a compatible manner with the CFL 
number in order to ensure monotonicity. We have omitted the formulas for the remaining cell 
boundaries, but they can be derived in a straightforward manner. Also note that even though the 
formulas are written along the boundary '/ -h 1/2, /, the use of the Gaussian quadrature formula 
will replace the axes'/ by Gauss points and therefore this subscript must be understood in that 
sense. It is important to notice that so far we are deriving expressions for the semi-discrete 
approximation of the system of conservation laws, however, the fluxes include the parameter 
At which arises from the averaging operators that originate these fluxes. Thus, in the spatial 
discretization of the system, the time stepping enters just as a parameter. At the end of the 
derivation of the scheme, when we present the time discretization of eq. (|4]), At will be considered 
as "marching parameter" in the sense that its inclusion in the formulas will generate an updated 
state in time. 

The function 1/2,7 " 0«+ 1/2,7(^^^1/2 j'^+1/2 P ^ limiter; a slight variation of the usual 
limiters has to be considered in this context since we use a centered flux instead of an upwind 
approach (the reader can refer to III Ch 13.] for more details); in our case we mainly use the 
SUPERB EE limiter, which on its centered version reads: 





2r 
1 

min {2, (pg -h (1 



if r < 0, 
ifO<r< ^, 
if I < r < 1, 



,)r\r>h 
5 



1-M 
1 + kl' 



(28) 



where c corresponds to the Courant number which depends on the problem. The Hmiter depends 
on the flow parameter r, which wiU be defined upon a physical quantity e of the system. Once e 
has been obtained from the discretized variables, left and right flow parameters are given by 



wo , -^f wo ^f,a/o , -^f 



^i+l/2J - R _ L ' ^+l/2j - P _L ' ^^^) 

1/2,7 ^^-+1/2,7 1/2,7 ^i+U2,j 



and finally, 



(A/+i/2,7 = min((A(rf^i/2j), ^(rf+i/2,7))- (30) 

The above described procedure starts with a set of averaged values and ends with a numerical 
approximation of the space operators involved in eq. ([2]). The resulting scheme is still continuous 
in time, and we conclude this section by discretizing this operator in a manner that is consistent 
with the choices that we have made in the generation of the space discretization operator. At a 
given starting time we begin by considering the semi-discrete scheme 

dQi j(t) 

^^=hj(Q), (31) 
at 

bringing the system to a final state f'^^ with a time stepping At. In order to preserve high- 
order and non-oscillatory properties in time, we consider the well-known family of explicit TVD 
Runge-Kutta schemes |8 |, in particular its third order version 

e;j^ = e«^. + ArL,,(Q«p, 02) 



^7 Q ^l,J Q ^l,J Q ^^J^^IJ ^ 



(33) 
(34) 



(35) 

We end this section with the inclusion of the source term. The source term appearing in eq. 
([3]), in the simplest case will not depend on space nor space derivatives, and therefore it can be 
averaged in space and solved in the same manner as the above presented time discretization, 
by replacing L/j(g^p by S(Q'lj). If we denote by the 2(At) the fully discrete operator that 
brings the system of conservation lawsd2]) At units ahead in time, and by (5 (AO the fully discrete 
operator that updates the source term ([Sjiin At units, we preserve, at least, second order accuracy 
in time by implementing a Strang splitting 1 14] in the form 

Q^y = G{At/2)2(At)6(At/2)0j. (36) 

If the source term does depends either on space or space derivatives, such a viscosity for instance, 
the averaging procedure will require the evaluation of the source term integral 

^iJ = irT- S{Q{x,z,t))dzdx. (37) 
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Proceeding in the same way as we did for the fluxes, we approximate this integral by a suitable 
double Gaussian quadrature, 

XGp ZGp 

where we make use of the same reconstruction procedure previously described in order to obtain 
values of Q inside every cell. 



3. Advective tests: Linear advection with space-dependent coefficients 

In this section we implement three diff'erent test cases based on the equation 

dtQ+d, (a(x, z)Q) + d, (b(x, z)Q) = 0, (39) 

which describes the evolution of a scalar quantity Q through a 2D domain fl, holding suitable 
initial and boundary conditions. The tests for this equation aims to recover the theoretically ex- 
pected second-order for the convergence rate of the scheme, to study variations on the weighting 
parameter oj in the flux and efl'ects of the limiter choice for sharp initial conditions. Test settings 
varies from one test to the other; though, one common aspect is the computation of the time 
stepping At, which once mesh parameters. Ax and Az, and the background flow (a, b) have been 
specified, is computed via 

At = CFL min( ^\ \, (40) 

\maxQ \a\ max^^ \b\i 

with CFL the classical Courant number which by default is set to CFL = 0.45, consistently 
chosen for a value of a; = 0.5 (there is a direct relation between the weighting parameter and the 
CFL number in order to preserve monotonicity of the scheme; see |[T3ll for precise details). A 
last common aspect the choice of flow parameter r, which is computed by taking ^ = g in eqs. 



3.L 2D Linear advection with constant coefficients 

The first case that we consider is a 2D linear advection equation with constant coefficients. 
We seta = \,b = 1,Q = [0, 1]^, and an initial profile given by 

2(x, z, 0) = sin(27rjc) sin(27rz), (41) 

together with periodic boundary conditions. We advect the initial profile for 10 periods, i.e., 
final time of the simulation is r = 10 [s]. Results for the convergence rates are shown in table [TJ 
expected second-order for smooth solutions is achieved in Li and Loo discrete norms, while an 
additional order is obtained for the Li norm, which is given by the fact that we are using a third 
order reconstruction in space. 
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Table 1: Convergence rates for the 2D linear advection problem with constant coefficients after t = I0[s] (10 periods). 



N 


Loo error 


Loo order 


Li error 


Li order 


50 


1.2637e-002 




2.005e-002 




100 


2.4316e-003 


2.4 


2.591 le-003 


3.0 


200 


6.1038e-004 


2.0 


3.5190e-004 


2.9 


400 


1.4790e-004 


2.0 


4.9012e-005 


2.9 



3.2. Swirling flow 

The second test that we present is a more stringent case, proposed in (15] . We make use of 
the setting proposed in 1 16 |: the velocity field is given by 

a = sin^(;rx) sin(27rz)g(t), b = - sin^(;rz) sin(2nx)g(t), (42) 

g(t) = cosintlT). (43) 
Again, is set to be the unit square, and the initial condition is taken as 

g(x,z,0) = i(l + cos(7rr), r = min(l, 4 V(^ - 0.25)^ + (z - 0.25)2). (44) 

Note that the velocity field vanishes at the boundary of the domain. Setting final time T = 5, we 
expect a maximal flow deformation at r/2 while, ideally, Q(x, z, T) = Q(x, z, 0). Figure [T] shows 
both exact initial (final) profile and maximal flow deformation at T/2 with Ax = Az = 0.005; the 
scheme preserves positivity of the initial profile at any time during simulation. Again, table |2] 
shows that second order of accuracy is reached as expected, although the convergence requires a 
larger number of elements; this is most likely due to the highly deformational nature of the flow. 



Table 2: Convergence rates for the swirling flow problem att = 5 [s]. 



N 


Loo error 


Leo order 


Li error 


Li order 


50 


7.1093e-001 




1.0719e-000 




100 


5.2363e-001 


0.4 


6.5837e-001 


0.7 


200 


2.4912e-001 


1.1 


2.2685e-001 


1.5 


400 


4.5618e-002 


2.4 


3.5510e-002 


2.7 



In this test we also include a study with low and high resolution simulations (Ax = 0.01 
and Ax = 0.005 respectively), where we try difl'erent values for the weighting parameter in the 
GFORCE flux, oj. Although the GFORCE flux used as low-order term in the FLIC approach 
is a first-order flux, (excepting for the case oj = I which corresponds to a Lax-Wendrofl' flux), 
increasing oj yields to a flux that has a behavior similar to a second-order flux; this is observed 
in figures [2] and [3j where noticeable increase on the accuracy of the final state is detected when 
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Figure 1: Swirling flow test problem. Initial condition (left) and maximal flow deformation a.tt = 2.5 [s] with 200 x 200 
elements (right). 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

X X 



switching to larger values of oj . Note that in every case the monotonicity of the solution is pre- 
served. There is an additional cost related with the increase of accuracy and preserving mono- 
tonicity at the same time, which is a decrease on the time stepping. For instance, according to 
1131, if OJ is between 0.5 and 1, 

At maxQ \a\ At max^ \b\ 
dx ' dz 

which implies that, for instance, when oj = 0.75, then CFL < 0.17. 

3.3. D osw ell fronto gene sis 

We conclude the advection tests with a kinematic frontogenesis problem, originally presented 
in 1 17 |. It is a standard test in atmospheric modelling, and allows us to assess the performance 
of the scheme in the treatment of sharp fronts; numerical experiments with this test case in the 
context of this article can be found in |[T8l[T9l. For this test, we set the domain = [-5, 5]^, and 
the flow is given by 

1 2 
a = -zf(r), b = xf(r), f(r) = -v(r), v(r) = v sech (r)tanh(r), (46) 

r 

r = ^|x^+z^, V = 2.59807. (47) 
Initial condition for this test is given by 



-l-\-oj 



2(JL> 



(45) 



Q(x,z,0) = tanh(^^y (48) 
9 



Figure 2: Swirling flow test problem. Low-resolution experiments with 100 x 100 elements att = 5 [s]. Varying values 
of CO in the GFORCE flux: From left to right, from top to bottom: oj = 0.25, 0.5, 0.75, 0.9. 
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Figure 3: Swirling flow test problem. High-resolution experiments with 200 x 200 elements att = 5 [s]. Varying values 
of CO in the GFORCE flux: From left to right, from top to bottom: co = 0.25, 0.5, 0.75, 0.9. 
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generating the following exact solution 



Q(x, z. t) = tank ^^-^ — . (49) 

The parameter 6 is related with the thickness of the front zone. We first present numerical exper- 
iments with a value of ^ = 1 , in order to generate a smooth solution and verify the convergence 
rates; table |3] shows that consistent convergence rates are obtained with the proposed scheme in 
both Loo and Li norms. For this case, final state profiles can be observed in figure [4j A second 
is study is performed with sharp fronts, taking 6 = 10"^. Figure [s] illustrates the capapcity of the 
scheme in tracking sharp fronts; it can be observed that no spurious oscillations are generated, 
despite the sharpness of the solution. Moreover, the scheme proves to be very robust with respect 
to the choice of the limiter, which is interesting as a common situation in the choice of limiters 
is that it highly depends on the problem. 



Table 3: Convergence rates for the Doswell frontogenesis problem at f = 4 [s], with S = I. 



N 



Loo error Loo order Li error Li order 



50 3.4786e-001 1.2719e-002 

100 1.1140e-001 1.6 3.5136e-003 1.9 

200 3.3302e-002 1.8 7.3045e-004 2.3 

400 5.47780e-003 2.6 1.0541e-004 2.8 



4. Convective tests: the Euler equations 

Throughout this section we consider a single model in order to study convective phenomena. 
Our starting point corresponds to the set of equations describing the evolution in time of a 2D 
dry air atmosphere |[20l O [Q. Imposing conservation of mass, momentum and energy, and 
considering eff'ects of gravity, together with neglecting friction and rotation eff'ects, leads to a set 
of 2D inviscid primitive equations for the atmosphere written in conservative form 



(50) 



where 



Q 



p 

pu 
pw 
p6 



r 



pu 
pu^ + 



puw 
pu6 



pw 
pwu 
pw 
pwO 



,2+ ^ 






-pg 





(51) 



In this system of equations p is the density of the fluid, u is the velocity in the x-direction, w is 
the velocity in the z-direction, P is the pressure and 6 is the potential temperature, which relates 
to the usual thermodynamic temperature T via 



= T 



r_ 
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-Rdlcp 



(52) 



Figure 4: Doswell frontogenesis problem. Results at t = 4[s] for smooth initial condition with S = I and 200 x 200 
elements. In the right bottom, a cut at x = 0. 




Figure 5: Doswell frontogenesis problem. Results at t = 4[s] for sharp initial condition with S = 10 ^ and 200 x 200 
elements. In the right bottom, a cut at x = with different limiters. 




The system is closed by the equation of state for an ideal gas 



P - C.m\ Co = (53) 

^0 

Model parameters are: the gravitational acceleration g = 9.81 [m^-"^], the atmospheric pressure 
at sea level Po = 10^ [Pa], the gas constant for dry air = 2^7[JK~^kg~^], the specific heat 
of dry air at constant pressure and volume Cp = 1004[JK~^kg~^], the specific heat of dry air at 
constant volume Cy = 717[JK~^kg~^] and its ratio y = CpjCy = 1.4. Additionally, by defining 
the Exner pressure n 



R/Cp 

(54) 



r_ 

To 

the expression for the total energy of the system (intemal+kinetic+potential) is given by 



e = Eint + EKin + Epot = CyOn + ^(^^ + w^) + gz, (55) 

which will be used as flow parameter in the limiter computation. For all the simulations in this 
section, the parameter co is set to 0.5, and once Ax and Az (although always Ax = Az) have been 
specified, the time stepping is selected according to 



At = CFL min 



Ax Az 



max^-v msixs. 



(56) 



with CFL = 0.4 and where Sx and are maximum characteristic speeds in the x and z direction 
respectively, 

Sx - max(w -\- Cs,u- c^), - max(w + c^, w - c^), (57) 

where Cs = yJdpP is the speed of sound in the fluid. In this section we study four test cases 
for the above presented set of equations. They all consist of hydrostatically balanced initial 
conditions plus a potential temperature perturbation. The first three cases are initialized with 
reference states for a neutral atmosphere, i.e., 

6 = 6-\-6\ 6 = 300[Kl (58) 
and the density p is initialized via the hydrostatic balance equation 

CpO^ = -g, (59) 

together with the ideal gas law 

P=^/r^^, (60) 

both evaluated at the reference state 6 = 6. The potential temperature perturbation, together with 
the initial velocity field and boundary conditions are specified for each problem. 
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4.1. Convective bubble in a neutral atmosphere 

This first test case, that has been previously addressed in ||2T]|22l|23l (among others), studies 
the behavior of a hot temperature bubble placed in a hydrostatically balanced neutral atmosphere 
at rest. As the perturbation is warmer than the background state, a buoyancy force will push it 
upwards and, as it starts rising, because of the same buoyancy eff'ect, it will start experiencing a 
deformation that will eventually develop into a mushroom-type of cloud. We use the scale and 
settings used in (H. The domain is D = [-10000, 10000] x [0, 10000], with Ajc = Az = 125[m] 
and the potential temperature perturbation is given by 

^, f2cos(f) L<1, L=-i-Vx2 + (z- 2000)2. (61) 
10 i.o.c. zUUU 

Simulation time has been set to 1000 [s], allowing the bubble to rise without hitting the top 
boundary; reflecting solid wall boundary conditions have been considered around the whole 
domain. 

It can be seen in figure[6]that the proposed scheme manages to reproduce the correct physical 
solution of the test; the bubbles rises from its initial state experiencing deformation that finally 
generates a mushroom cloud. Note that symmetry is preserved during the simulation, no spurious 
oscillations are detected and the final front location is essentially the same as shown in ll24ll . also 
matching in the vertical velocity, as shown in figure |7] It has to be remarked that most of the 
numerical models used to approximate such experiment need to include a damping eff'ect in 
order to obtain a grid converged solution, which is not our case. We also assess the performance 
of the scheme in this test by the energy plot in figure [Sl it illustrates the quantities present in 



eq. ( [55] ) (we consider the the quantity pe rather than e), and it can be seen the rise of the kinetic 
energy of the system initially at rest together with a decrease of the internal and potential energy 
contributions to the system, all this with a constant total energy, as it is expected for a closed 
system without physical dissipation mechanisms. 

4.2. Interaction between hot and cold bubbles in a neutral atmosphere 

We turn our attention to a variation of the previously presented case. Similar tests have being 
previously performed in ll2T1l : we have modified it in order to consider the same scales as in 
our first case, but also to include some kinetic eff'ects to test symmetry. The perturbation in the 
potential temperature contains a hot but also a cold bubble; the expressions are given by 

e = e + e[+e2, (62) 

^ /lOcos(f) Li<l,^ ^^^|_i5cos(f ) L,<1,^ ^^^^ 
i.o.c. ' |0 i.o.c. 



Li = ^Vx2 + (^_ 2000)2, L2= ^V-^2 + (^ _ 8000)2. (64) 

(65) 



We also expect in this case that the hot bubble will rise, but on the other hand the cold bubble 
should fall, and as they are placed along the same vertical axis, they will collide and interact 
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Figure 6: Convective bubble in a neutral atmosphere test case; potential temperature colormaps. From left to right, top 
to bottom: initial state and Results at t = 300, 600 and 1000 [s], with Ax = Az = 125[m], 160 x 80 elements. 
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Figure 7: Convective bubble in a neutral atmosphere test case; velocity field colormaps at f = 600 and t = 1000 [s], with 
Ax = Az = 125[m], 160 x 80 elements. Left: horizontal velocity u. Right: vertical velocity w. 
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creating eddy patterns. We consider solid wall boundary conditions at the bottom and the top 
of the domain, and periodic boundary conditions at the lateral extremes, as we also include an 
initial horizontal velocity u = 20[ms~^] to test the capacity of the scheme to preserve symmetries 
in the presence of horizontal translation. Figure [9] illustrates the initial setting, figures 10 and 



11 exhibiting the evolution of the system with a resolution of Ax = Az = 125[m] up to 1000[5'] 
The results reflect the proper physical solution, the rise of the hot bubble together with the fall 
of the cold bubble until the collision is produced, and the consequent eddy generation due to 
the interaction of the perturbations. The rotational behavior can be also noticed in the velocity 
plots, where eddies are also formed. The final state at 1000[^] is symmetric with respect to the 
axis X = 0, which is coherent with the initial horizontal velocity. No spurious oscillations are 
observed. The energy plot in figure [12] shows an increment in the kinetic energy of the system, 
interaction between internal and potential energy, but shows almost no artificial diff'usion for the 
total energy of the system. 

4.3. Density current 

The third case is a popular test case in atmospheric modelling (see |f25l[23l|26l[27l); in the 
domain VI = [0, 20000] x [0, 6000], with a system initially at rest, a cold bubble perturbation is 
added to the reference state, 

^^j-7.5(C0S(,Z.,.,) L,U^ ,,,„^|,|^,, 

With solid boundary conditions, the cold bubble drops until it hits the boundary, generating hor- 
izontal shear displacement, end eventually developing Kevin-Helmholtz rotors. Final simulation 
time is t = 900[^]. It has been previously reported in ||25]|23l that most numerical methods 
require to take into account viscosity eff'ects in order to generate a grid-convergent solution; nev- 
ertheless, there exists methods that manage to accurately simulate this test case with a inviscid 
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Figure 9: Colormap of the initial potential temperature for the hot and cold bubbles test case. 
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Figure 10: Hot and cold bubbles test case.Results at t = 180 and 250 [s], with Ax = Az = 125[m], 160 x 80 elements. 
Left: colormap of the potential temperature. Right: vector plot of the velocity field. 
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Figure 11: Hot and cold bubbles test case.Results at t = 500 and 1000[s], with Ax = Az = 125[m], 160 x 80 elements. 
Left: colormap of the potential temperature. Right: vector plot of the velocity field. 




Figure 12: Normalized energy (with respect to the initial total value) of the system for the hot and cold bubbles test case. 
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set of equations f28l . When a viscous stress term is included in the scheme, a reconstruction step 
has to be included in the source term approximation and second spatial derivatives are computed 
from the reconstructed polynomials. In our case, we include simulations for the inviscid model, 
but we also add a viscous source term to the r.h.s. of the system ([5Q|), with 





pK(dlu + d^u) 
pK(dlw + d^w) 
pK(dl6 + d^e) 



(67) 



We first study the eff'ect of resolution in both inviscid and viscous simulations. In figures 
[13] and [14] it can be seen final time results for variable resolution. In both cases, there is a clear 
convergent behavior of the solution, and increasing resolution provides a better insight into the 
development of Kevin-Helmhotlz rotors; at the highest tested resolution, with Ax = Az = 50[m], 
it can be clearly appreciated the generation of three rotors, which is similar to the results obtained 
in the aforementioned references for this test case for viscous simulations (including the diff'used 
aspect of the eddy nearest to the front). Table [4] shows the extreme values obtained for both 
simulations at high-resolution, which are in accordance to the range of values previously obtained 
by other authors; in particular, the front location, which is a relevant quantity in this test case, 
coincides the results published in ||23]|25]|271. The energy plots shown in figure [T5| and [T6| are 
qualitatively comparable to the one presented in 1261 (even though is not exactly the same density 
current test case), exhibiting a sustained increment in the kinetic energy of the system, while both 
internal and potential energy decrease throughout the simulation. There are some diff'erences 
though, when it comes to analyze the conservation of the total energy: the inviscid system is 
completely conservative, while the inclusion of diff'usion alters this property and decreases in 
time. 



Table 4: Extreme values for the density current test case with Ax = Az = 50[m] at t = 900[^]. 



Model 


6' . 

mm 


max 


^min 


Umax 


^min 


^max 


Front location 


Inviscid 
Viscous 


-11.7296 
-13.1851 


0.8950 
1.0209 


-18.6235 
-17.2731 


32.2655 
30.5931 


-18.9703 
-15.5440 


22.6452 
15.7718 


1.498x104 
1.503x10^ 



4.4. Convective bubble in a stable atmosphere 

We conclude our study with a test that has been previously presented in ||22| . In the domain 
VI = [0, 40000] X [0, 15000], we consider a stable atmosphere with an initial potential temper- 
ature vertical gradient of 4[Kkm~^] with a mean ground level value of 300[^]. After vertical 
integration of the gradient, a positively stratified reference state for the potential temperature is 
obtained, and density is initialized via hydrostatic balance. We add a warm potential temperature 
perturbation bubble of the form 

6^ = 1 ^=^V^' + fe- 2750)2. (68) 

10 i.o.c. 25UU 

Simulation time is set to ^ = 600[^], we use solid wall boundary conditions at the top and the 
bottom of the domain, and open boundary conditions for the lateral extremes. We execute model 
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Figure 13: Density current test case. Results after 15 min. From top to bottom: potential temperature colormap with 
Ax = Az = 200, 100and50[m]. Left: experiments without viscosity. Right: results with Fickian viscosity with parameter 
y = 15[m^ s-^]. 
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Figure 14: Density current test case. Results after 15 min and Ax = Az = 50[m], 400 x 120 elements. Top: horizontal 
velocity. Bottom: vertical velocity. Left: experiments without viscosity. Right: results with Fickian viscosity with 
parameter v = 75 [m^ s~^]. 




Figure 16: Normalized energy (with respect to the initial total value) of the system for the density current test case with 
viscosity. 
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runs with Ax = Az = 500[m] and Ax = Az = 250[m]. Figure 17 shows high and low-resolution 
results at final time of simulation. Low-resolution results are in accordance to what is presented 
in 1 29 1; since the temperature increases with altitude, a mitigation of the buoyancy is expected 
while horizontal spreading of the perturbation occurs. High-resolution experiments exhibits the 
same behavior, although there is a variation in the extremal values as it can be seen in table [5] In 
both cases symmetry with respect to x = is preserved. Regarding energy conservation, figure 
18 shows that total energy is preserved while potential, internal and kinetic energy oscillate with 
respect the initial state. 



Table 5: Extreme values for the convective bubble in a stable atmosphere test case at r = 600[^]. 



Resolution 


& . 

mm 


max 


^min 


Umax 


^min 


^max 


Ax = 500 [m] 
Ax = 250 [m] 


0.0748 
-0.0251 


10.3438 
10.3890 


-1.9717 
-3.0724 


1.9717 
3.0685 


-1.6539 
-2.3212 


1.4501 
2.1149 



5. Summary and outlook 

We have presented a second-order, non-oscillatory scheme for the resolution of relevant ad- 
vective and convective atmospheric phenomena. The method makes use of a WENO recon- 
struction procedure for accurate extrapolation of boundary and inner cell values, together with 
a centered-limited approach for the flux calculation. Even though, for this class of problems, 
schemes based on up winding considerations are usually preferred over centered approaches, the 
proposed scheme performs well in an extensive set of test problems. The theoretically expected 
second order was reached whenever an analytic solution was available, and in its absence the 
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Figure 17: Convective bubble in a stable atmosphere. Results after 10 min. From top to bottom: potential temperature, 
horizontal velocity and vertical velocity. Left: low resolution results with Ax = Az = 500[m]. Right: high resolution 
results with Ax = Az = 250 [m]. 
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Figure 18: Normalized energy (with respect to the initial total value) of the system for convective bubble in a stable 
atmopshere. 
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scheme proved to be grid-convergent with a performance similar to the currently available algo- 
rithms. The accuracy of the method can be numerically increased by switching the flux parameter 
oj, with a cost associated to the time stepping in order to preserve monotonicity. The scheme also 
showed robustness with respect to the limiter choice in a sharp front, which also illustrates the ro- 
bustness of the reconstruction procedure. In the convective experiments, the scheme managed to 
reproduce the correct physical behavior, while tracking fronts in the correct position and without 
generating spurious oscillations. It has also showed robustness with respect of the inclusion of 
viscosity: unless it is physically relevant, the current evidence seems to indicate that the method 
is able to perform in a consistent way without the need of viscosity, which is a consequence of 
the continuous enforcement of the non-oscillatory character of the scheme, present in both time 
and spatial discretizations. 

We point some open issues to address. The removal of the splitting approach, by trying 
to incorporate source terms efl'ects into the flux calculation could make the method far more 
efliicient. In the splitting, context, other alternatives for viscosity calculation should be explored, 
as in the current scheme additional reconstruction steps (one of the most expensive parts of the 
code) are required. Nevertheless, we conclude that the studied scheme possess a great level of 
applicability in the atmospheric modelling framework; further enhancements on its formulation 
and performance could increase its potential as a robust, high-resolution method. 
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